# TODO: create PM image from aod and temp (PM regression)
# 
# Author: phamha
###############################################################################

#import library
library(gstat)
library(base)
library(RPostgreSQL)
library(stringr)
library(raster)
library(gdalUtils)
library(rgdal)
library(rPython)


#host_name = '172.16.81.252'
host_name = '112.137.129.222'
database_name = 'apom'
user_name = 'postgres'
password = 'postgres'

min_aod_mod = 0.0015
max_aod_mod = 1.165031427
min_temp_mod = 291.8899935
max_temp_mod = 317.5399929

min_aod_myd = 0.057872928
max_aod_myd = 1.509154063
min_temp_myd = 281.2699937
max_temp_myd = 316.6299929

min_aod_npp = 0.065608211
max_aod_npp = 1.801370621
min_temp_npp = 281.0396729
max_temp_npp = 326.4881287

aod_scale_mod = 0.00100000004749745
aod_scale_myd = 0.00100000004749745
aod_scale_npp = 1


res_folder_mod = "res/SatResampMOD04"
prod_folder_mod = "prod/ProdMODPM"

res_folder_myd = "res/SatResampMYD04"
prod_folder_myd = "prod/ProdMYDPM"

res_folder_npp = "res/SatResampleViirsAOT"
prod_folder_npp = "prod/ProdVIIRS"

met_folder_mod = "/var/www/html/MODIS/HaPV/MetTiff/"
met_folder_npp = "/var/www/html/MODIS/HaPV/MetNPP/"
shape_file = "/var/www/html/MODIS/HaPV/BDNEN/VNM_adm0.shp"
tif2raster_file = "/var/www/html/MODIS/HaPV/tif2rasref.py"
create_aqi_file = "/var/www/html/MODIS/HaPV/MODIS_CLIP_PM_AQI.py"
create_png_file = "/var/www/html/MODIS/HaPV/create_png.py"
create_aot_file = "/var/www/html/MODIS/HaPV/aot_processing.py"




insert_all_query = "INSERT INTO %s(aqstime, rasref, filename, filepath, gridid,pmid,max,min,avg, type,sourceid) VALUES ('"
get_time_query = "select aqstime from %s where filename like "
listfile_query = "SELECT mod04.aqstime, mod04.filename,mod04.filepath,mod07.filename as temp_filename,mod07.filepath as temp_filepath FROM res.satresampmod04 as mod04 inner join res.satresampmod07temperature as mod07 ON (mod04.aqstime = mod07.aqstime) where mod04.aqstime > '2014-07-25 00:00:00'::timestamp order by mod04.aqstime"


regress_predict = function(sate_data,model_type,cook_type,aod,temp,avg_temp,avg_rh,avg_preci){
	if(sate_data=="mod"&&model_type=="linear"&&cook_type=="4np"){
		pm25 = 22.58132724* aod + 0.08452807 * temp + (-2.98371028)*avg_temp + (-0.58041657)*avg_rh + 0.64525692 * avg_preci + 113.98273058
	}
	if(sate_data=="mod"&&model_type=="linear"&&cook_type=="fdist"){
		pm25 = 29.3009026972483* aod - 0.0121947942915479 * temp + (-4.05178592218209)*avg_temp + (-1.17578020320201)*avg_rh + 1.01861037054962 * avg_preci + 213.045850613292
	}
	if(sate_data=="myd"&&model_type=="linear"&&cook_type=="4np"){
		pm25 = 25.6965586 * aod + 0.1525928 * temp + (-2.1894375)*avg_temp + (-0.6837351)*avg_rh + 0.4189560 * avg_preci + 70.1740458
	}
	if(sate_data=="myd"&&model_type=="linear"&&cook_type=="fdist"){
		pm25 = 30.5744605532908 * aod + 0.118323308557864 * temp + (-2.47362345170908)*avg_temp + (-0.795065245241148)*avg_rh + 0.45323741215544 * avg_preci + 94.9758204191974
	}
	if(sate_data=="npp"&&model_type=="linear"&&cook_type=="4np"){
		pm25 = 14.9045673 * aod + 0.0647682902577907 * temp + (-2.91326445481888)*avg_temp + (-1.20281335981164)*avg_rh + 0.446332747336052 * avg_preci + 154.899079958179
	}
	
	return(pm25)
}
#get data from DB
getDataFromDB = function(sql_command){
	driver = dbDriver("PostgreSQL")
	connect = dbConnect(driver,dbname = database_name,host = host_name,port=5432,user = user_name,password= password)
	rs = dbSendQuery(connect,sql_command)
	data=fetch(rs,n=-1)
	dbDisconnect(connect)
	dbUnloadDriver(driver)
	return (data)
}
#insert data to DB
insertDataToDB = function(sql_command){
	driver = dbDriver("PostgreSQL")
	connect = dbConnect(driver,dbname = database_name,host = host_name,port=5432,user = user_name,password= password)
	rs = dbSendQuery(connect,sql_command)
	dbDisconnect(connect)
	dbUnloadDriver(driver)
}

createKrigingImage = function(regressPm_file){
	
	#PM values
	pmRaster=raster(regressPm_file)
	pm=values(pmRaster)
	corxy=coordinates(pmRaster)
	x=corxy[,'x']
	y=corxy[,'y']
	
	
	totalCell=length(pmRaster)
	cell = c(1:totalCell)
	
	table=data.frame(cell,x,y,pm)
	newTable=table
	testTable=subset(table,pm<0)
	trainTable=subset(table,pm>=0)
	
	#caculate variogram
	empiVario=variogram(pm~1,locations=~x+y,data=trainTable)
	
	#sph fit
	sill=min(empiVario$gamma)
	#sphModel=vgm(psill=sill,model="Sph",nugget=0,range=min(empiVario$dist))
	sphModel=vgm(model="Sph",nugget=0,range=1)		
	sphFit=fit.variogram(empiVario,sphModel)
	
	universal_result=krige(id="pm",formula=pm~x+y,data=trainTable,newdata=newTable,model=sphFit,locations=~x+y)
	
	#edit tiff
	universalPMRaster=pmRaster
	universalPMValue=universal_result[,3]
	universalPMRaster[1:totalCell]=universalPMValue
	
	
	#edit error tiff
	errorPMRaster=pmRaster
	errorPMValue=universal_result[,4]
	errorPMRaster[1:totalCell]=errorPMValue
	
	#save uk result to tiff
	uk_file = str_replace(regressPm_file,"rg.tif","uk.tif")
	writeRaster(universalPMRaster,filename=uk_file,format="GTiff",overwrite=TRUE)
	gdal_rasterize(shape_file,uk_file,b=1,i=TRUE,burn=-9999,l="VNM_adm0")
	
	#set n/a value
	new_uk_raster = raster(uk_file)
	new_uk_value = values(new_uk_raster)
	new_uk_value[new_uk_value==-9999]<-NA
	new_uk_value[new_uk_value<0]<-0
	new_uk_raster[1:totalCell] = new_uk_value

	writeRaster(new_uk_raster,filename=uk_file,format="GTiff",overwrite=TRUE)
	
	
	#save uk error to tiff
	error_file = str_replace(regressPm_file,"rg.tif","error.tif")
	writeRaster(errorPMRaster,filename=error_file,format="GTiff",overwrite=TRUE)
	gdal_rasterize(shape_file,error_file,b=1,i=TRUE,burn=-9999,l="VNM_adm0")
	print(uk_file)
	
}
create_pm_image = function(sate_data,model_type,cook_type,aod_file,temp_file,source_id){
	#print(aod_file)
	#print(temp_file)
	if(sate_data=="mod"){
		type = 0
		insert_query = sprintf(insert_all_query, "prodpm.prodmodispm_vn_collection0")
		met_folder = met_folder_mod
		aod_scale = aod_scale_mod
		time_query = sprintf(get_time_query, "res.satresampmod04")
		res_folder = res_folder_mod
		prod_folder = prod_folder_mod
		min_aod = min_aod_mod
		max_aod = max_aod_mod
		min_temp = min_temp_mod
		max_temp = max_temp_mod
		
	}
	if(sate_data=="myd"){
		type = 1
		insert_query = sprintf(insert_all_query, "prodpm.prodmodispm_vn_collection0")

		met_folder = met_folder_mod
		aod_scale = aod_scale_myd
		time_query = sprintf(get_time_query, "res.satresampmod04")
		res_folder = res_folder_myd
		prod_folder = prod_folder_myd
		min_aod = min_aod_myd
		max_aod = max_aod_myd
		min_temp = min_temp_myd
		max_temp = max_temp_myd
	}
	if(sate_data=="npp"){
		type = 2
		insert_query = sprintf(insert_all_query, "prodpm.prodviirspm_vn_collection0")
		met_folder = met_folder_npp
		aod_scale = aod_scale_npp
		time_query = sprintf(get_time_query, "res.satresampviirs")
		min_aod = min_aod_npp
		max_aod = max_aod_npp
		min_temp = min_temp_npp
		max_temp = max_temp_npp
		res_folder = res_folder_npp
		prod_folder = prod_folder_npp
	}
	filename = basename(aod_file)
	filename = str_replace(filename,".tif","")
	
	#start_index = regexpr("M[^M]*$",aod_file)
	#end_index = regexpr(".tif",aod_file) - 1
	#file_name = substr(aod_file,start_index,end_index)
	
	time_query = paste(time_query,file_name,"%'",sep="")
	data = getDataFromDB(time_query)
	
	mod04_aqstime = data$aqstime[1]
	aqstime = strptime(mod04_aqstime,format="%Y-%m-%d %H:%M:%S")
	aqstime = aqstime + 25200	
	month = format.Date(aqstime,"%m")
	year = format.Date(aqstime,"%Y")
	
	#mask VietNam base on shapfile
	aod_mask_file = str_replace(aod_file,".tif","_mask.tif")
	file.copy(aod_file,aod_mask_file)
	
	gdal_rasterize(shape_file,aod_mask_file,b=1,i=TRUE,burn=-9999,l="VNM_adm0")
	aod_dataset = raster(aod_mask_file)
	aod_dataset[aod_dataset[] == -9999] <- NA
	writeRaster(aod_dataset,filename=aod_mask_file,format="GTiff",overwrite=TRUE)
	
	aod_data = values(aod_dataset)
	aod_data = aod_data * aod_scale
		
	corxy = coordinates(aod_dataset)
	x = corxy[,'x']
	y = corxy[,'y']
	
	
	temp_mask_file = str_replace(temp_file,".tif","_mask.tif")
	file.copy(temp_file,temp_mask_file)
	gdal_rasterize(shape_file,temp_mask_file,b=1,i=TRUE,burn=-9999,l="VNM_adm0")
	temp_dataset = raster(temp_mask_file)
	temp_dataset[temp_dataset[] == -9999] <- NA
	writeRaster(temp_dataset,filename=temp_mask_file,format="GTiff",overwrite=TRUE)
	
	temp_data = values(temp_dataset)
	if(sate_data!="npp"){
		temp_data = (temp_data + 15000) * 0.00999999977648258
	}
	
	avg_temp_file  = paste(met_folder,"temp",as.numeric(month),".tif",sep="")
	avg_rh_file    = paste(met_folder,"rh",as.numeric(month),".tif",sep="")
	avg_preci_file = paste(met_folder,"preci",as.numeric(month),".tif",sep="")
		
		
	avg_temp_dataset  =  raster(avg_temp_file)
	avg_rh_dataset    =  raster(avg_rh_file)
	avg_preci_dataset =  raster(avg_preci_file)
		
	avg_temp_data  = values(avg_temp_dataset) 
	avg_rh_data    = values(avg_rh_dataset) 
	avg_preci_data = values(avg_preci_dataset)
		

	pm25 = regress_predict(sate_data,model_type,cook_type,aod_data,temp_data,avg_temp_data,avg_rh_data,avg_preci_data)
	
	total_pixel = sum(!is.na(aod_data))
	if(sate_data!="npp"){
		ratio = total_pixel/2024*100
	}else{
		ratio = total_pixel/11180*100
	}
	print(paste("Pixel:",total_pixel,"Cloud ratio:",ratio,sep=" "))
	
	if(ratio>=30){
		# create regression image
		table = data.frame(x,y,aod_data,temp_data,avg_temp_data,avg_rh_data,avg_preci_data,pm25)
		table$pm25[table$aod_data<min_aod|table$aod_data>max_aod|table$temp_data<min_temp|table$temp_data>max_temp]<-NA
		og_raster = aod_dataset
		totalCell = ncell(og_raster)
		og_raster[1:totalCell] = table$pm25
		
		if(sate_data!="npp"){
			pm_file = str_replace(aod_file,".hdf_DT_10km.tif","_rg.tif")
			pm_file = str_replace(aod_file,".hdf_CB_10km.tif","_rg.tif")
			
		}else{
			pm_file = str_replace(aod_file,".tif","_rg.tif")
			
		}
		pm_file = str_replace(pm_file,res_folder,prod_folder)
		
		print (pm_file)
	
		#mid_index = regexpr("M[^M]*$",pm_file)
		#out_path = substr(pm_file,0,mid_index-1)
		out_path = dirname(pm_file)
		dir.create(path=out_path,showWarnings = FALSE,recursive = TRUE,mod="777")
	
		writeRaster(og_raster,filename=pm_file,format="GTiff",overwrite=TRUE)
		#gdal_rasterize(shape_file,pm_file,b=1,i=TRUE,burn=-9999,l="VNM_adm0")
		print ("Create regression image good!");

		#create pm image
		createKrigingImage(pm_file)
		print ("Create Kriging image good!");
	
		uk_file = str_replace(pm_file,"rg.tif","uk.tif")
		python.assign("raster_file",uk_file)
		python.load(tif2raster_file)
		raster_ref = python.get("raster_ref")
		
		uk_raster = raster(uk_file)
		uk_value = values(uk_raster)
		uk_value = uk_value[uk_value!=-9999]
		max_value = max(uk_value,na.rm = TRUE)
		min_value = min(uk_value,na.rm = TRUE)
		avg_value = mean(uk_value,na.rm = TRUE)
	
		start_index = regexpr("apom/prod",uk_file)
		
		if(sate_data!="npp"){
			mid_index = regexpr("M[^M]*$",uk_file)
			
		}else{
			mid_index = regexpr("G[^G]*$",uk_file)
			
		}
		#mid_index = regexpr("M[^M]*$",uk_file)
		end_index = nchar(uk_file)	
		uk_file_path = substr(uk_file,start_index,mid_index-1)
		uk_file_name = basename(uk_file)
		#uk_file_name = substr(uk_file,mid_index,end_index)
	
		aqstime2 = aqstime - 25200
		query = paste(insert_query,aqstime2,"'::timestamp, '",raster_ref,"'::raster, '",uk_file_name,"', '",uk_file_path,"', 1, 1, ",max_value,", ",min_value,", ",avg_value,", ",type,", ",source_id,")",sep="")
		
		#print(query)
		#insert pm images to database
		#insertDataToDB(query)
		print ("Insert to database good!");
		
		#create png
		print ("Tao anh PNG");
		python.load(create_png_file)
		
		#python.call("CreatePM25PNGImg",basename(uk_file),dirname(uk_file))

		#create aqi
		print ("Den phan cua chuc");
		print ("Cat anh AQI");
		python.assign("pmname",uk_file_name)
		python.assign("sat_type",type) # sat_type is mod:0, myd:1
		#python.load(create_aqi_file)

		# cat anh AOT
		print ("Cat anh AOT");
		#python.assign("tifname", basename(commandArgs(TRUE)[2])) # load the MOD04/MYD04 filename
		#python.assign("sat_type",type) # sat_type is mod:0, myd:1	
		#python.assign("sourceid",commandArgs(TRUE)[4]) # sourceid
		
		#print(aod_file)
		#print(type)
		#print(source_id)
		#python.load(create_aot_file)
		#python.call("ExecuteAll", basename(aod_file), type, source_id)
		
		
		
	}
			
}


# Test example
sat_data = "npp"
source_id = 0

aod_file = "/apom_data/apom/res/SatResampleViirsAOT/2014/GAERO-VAOOO_npp_d20140124_t0606539_e0612325_b11618_c20140725113835527415_noaa_ops_550_resample/GAERO-VAOOO_npp_d20140124_t0606539_e0612325_b11618_c20140725113835527415_noaa_ops_550_resample.tif"
temp_file= "/apom_data/apom/res/SatResampVIIRSTemp/2014/GMTCO-VLSTO_npp_d20140124_t0606539_e0612325_b11618_c20150119071449593020_noaa_ops/GMTCO-VLSTO_npp_d20140124_t0606539_e0612325_b11618_c20150119071449593020_noaa_ops_LST_resample.tif"

# mod/myd/npp , aot file, temp file
# Test command
create_pm_image("npp","linear","4np",aod_file,temp_file, source_id)
print("done")


# Real command
# command 1: type (mod/myd)
# command 2: mod04
# command 3: mod07
# command 4: sourceid (0: NASA, 1: UET)
#create_pm_image(commandArgs(TRUE)[1],"linear","4np",commandArgs(TRUE)[2],commandArgs(TRUE)[3],commandArgs(TRUE)[4])



## listfile_query = "SELECT mod04.aqstime, mod04.filename,mod04.filepath,mod07.filename as temp_filename,mod07.filepath as temp_filepath FROM res.satresampmod04 as mod04 inner join res.satresampmod07temperature as mod07 ON (mod04.aqstime = mod07.aqstime) where mod04.aqstime > '2014-07-25 00:00:00'::timestamp and mod07.filename like '%10km%' order by mod04.aqstime"
## data = getDataFromDB(listfile_query)
## total_record = nrow(data)
## #print(data)

## for(i in 1:total_record){
## 
##     aod_filename = str_trim(data$filename[i])
##     aod_path = str_trim(data$filepath[i])
##     aod_file = paste(data_folder,aod_path,aod_filename,sep = "")
## 
## 
##     temp_filename = str_trim(data$temp_filename[i])	
##     temp_path = str_trim(data$temp_filepath[i])
##     temp_file = paste(data_folder,temp_path,temp_filename,sep = "")
##     create_pm_image(sate_data,model_type,cook_type,aod_file,temp_file,source_id)
##     create_pm_image("mod","linear","4np",aod_file,temp_file)
## 
## 
## }


